Logistic map and Lyapunov exponent

We will plot the bifurcation diagram of a nonlinear dynamical system and use the Lyapunov exponent to identify chaotic behavior. The logistic map is deterministic : the same parameter and initial condition define the same mathematical trajectory. Nevertheless, in the chaotic regime, tiny differences in the initial condition—including measurement uncertainty and floating-point rounding—grow exponentially. Long-term prediction then becomes practically impossible even though the governing equation is deterministic.

A chaotic dynamical system is highly sensitive to initial conditions; small perturbations at any given time yield completely different trajectories. The trajectories of a chaotic system tend to have complex and unpredictable behaviors.

Many real-world phenomena are chaotic, particularly those that involve nonlinear interactions among many agents (complex systems). Examples can be found in meteorology, economics, biology, and other disciplines.

In this recipe, we will simulate a famous chaotic system: the logistic map. This is an archetypal example of how chaos can arise from a very simple nonlinear equation. The logistic map models the evolution of a population, taking into account both reproduction and density-dependent mortality (starvation).

We will draw the system's bifurcation diagram, which shows the possible long-term behaviors (equilibria, fixed points, periodic orbits, and chaotic trajectories) as a function of the system's parameter. We will also compute an approximation of the system's Lyapunov exponent, characterizing the model's sensitivity to initial conditions.

The logistic map is defined by the recursive application of the logistic function:

$$x_{n+1}= r x_n (1-x_n)$$

It is a discrete demographic model, where $x_n$ is a number between zero and one representing the ratio of the existing population to the maximum possible population.

The population is expected to increase at a rate proportional to the current population when the population size is small. The proportionality is a parameter $r$. However, when population grows large it starts to overuse available resources, and starvation occurs, hence the growth rate will decrease at a rate proportional to the value obtained by taking the theoretical "carrying capacity" of the environment (here 1.0) less the current population.

The usual values of interest for the parameter $r$ are in the interval $[0,4]$, for which the map sends $[0,1]$ into itself. For $r>4$, some initial conditions in this interval produce values outside it, including unphysical negative population sizes.

For a bifurcation diagram, most generic initial conditions with $0<x_0<1$ approach the same long-term attractor for a given $r$. We choose $x_0=10^{-5}$. Exceptional initial conditions exist—for example, $x_0=0$ remains zero forever—and individual chaotic trajectories are extremely sensitive to the precise value of $x_0$.

Notice that logistic map could also be expressed as $f(x) = r x (1-x)$ and we are looking for $f(f(f(...f(x_0)...)))$.

We will also calculate the Lyapunov exponent, which characterizes the rate of separation of infinitesimally close trajectories.

Two trajectories in phase space with initial separation $\delta Z_0$ diverge at a rate given by $$|\delta Z(t)|\approx e^{\lambda t}|\delta Z_0|.$$ If $\lambda>0$, nearby trajectories separate exponentially and the system is chaotic. If $\lambda<0$, nearby trajectories converge toward a stable fixed point or periodic orbit. Values near $\lambda=0$ occur at transitions and bifurcations.

The formal definition of Lyapunov exponent is $$\lambda=\lim_{t\rightarrow \infty}\lim_{|\delta Z_0|\rightarrow 0}\frac{1}{t}\ln(\frac{|\delta Z(t)|}{|\delta Z_0|})$$

For discrete map $x_{n+1}=f(x_n)$

$$\delta x_{n+1} \equiv x^1_{n+1}-x^0_{n+1} = f(x^1_{n})-f(x^0_{n}) \approx f'(x_{n})\delta x_n \approx \delta x_0 \prod_{i=0}^n f'(x_i) $$

After taking the logarithm, and dividing by the number of steps, we get $$\lambda = \lim_{n\rightarrow\infty}\frac{1}{n}\ln(\frac{|\delta x_{n}|}{|\delta x_0|})=\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=0}^{n-1} \ln|f'(x_i)|$$

For logistic map, the Lyapunov exponent is hence: $$\lambda =\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=0}^{n-1} \ln|r (1-2 x_i)|$$

For the algorithm, we first perform Ntransient iterations that are discarded, allowing the trajectory to approach its long-term attractor. We then perform Nmeasure additional iterations, use them to calculate the Lyapunov exponent, and save the last Nsave points for the bifurcation diagram.

The parameter r will be discretized with linear mesh of Npoints between ext[0] and ext[1].

The trajectories corresponding to different values of r are independent. We therefore parallelize the outer loop over the parameter index with Numba's prange. We create an array r[0:Npoints], store the final orbit points in data[0:Nsave,0:Npoints], and calculate one Lyapunov exponent for each value of r.

Here is the code:

import numpy as np
from numba import njit, prange

@njit 
def logistic(r, x):
    return r * x * (1 - x)

@njit(parallel=True)
def GiveFewPoints(ext, Npoints, Ntransient, Nmeasure, Nsave):
    """ ext     : contains limit for value of r: [r_start,r_end]
        Npoints : is number of r's between ext[0] and ext[1]
        Ntransient: iterations discarded before measurement
        Nmeasure : iterations used to measure the long-term behavior
        Nsave   : the last few steps should be saved. Hence Nsave should be 
                  much smaller than Nmeasure
    """
    r = np.linspace(ext[0], ext[1], Npoints)
    data = np.zeros((Nsave, Npoints))
    lyapunov = np.zeros(Npoints)

    # Each value of r defines an independent trajectory.
    for ir in prange(Npoints):
        rr = r[ir]
        x = 1e-5

        # Approach the long-term attractor before making measurements.
        for _ in range(Ntransient):
            x = logistic(rr, x)

        total = 0.0
        for i in range(Nmeasure):
            # The derivative f'(x_i)=r(1-2x_i) is evaluated before
            # advancing to x_{i+1}, matching the mathematical formula.
            total += np.log(np.abs(rr * (1.0 - 2.0*x)))
            x = logistic(rr, x)

            if i >= Nmeasure - Nsave:
                data[i - (Nmeasure - Nsave), ir] = x

        lyapunov[ir] = total / Nmeasure

    return r, data, lyapunov

We sample Npoints=1000 parameter values. For each trajectory, we discard Ntransient=1000 iterations, measure the following Nmeasure=9000 iterations, and save the final Nsave=300 generations.

ext = np.array([1.0, 4.0])
r, data, lyapunov = GiveFewPoints(ext, Npoints=1000, Ntransient=1000, Nmeasure=9000, Nsave=300)

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

We will use two panels that share the same $r$ axis. The upper panel is the bifurcation diagram: for each value of $r$, it shows the last Nsave iterates after the transient has been discarded. The lower panel shows the corresponding Lyapunov exponent.

The format string 'k,' draws very small black pixels, and alpha=0.2 makes overlapping points easier to see. The horizontal line $\lambda=0$ separates contracting behavior ($\lambda<0$) from sensitive, chaotic behavior ($\lambda>0$).

import matplotlib.pyplot as plt
%matplotlib inline

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(8, 9))

# Plot every saved iterate in one call.
ax1.plot(np.tile(r, data.shape[0]), data.ravel(), 'k,', alpha=0.2)
ax1.set_ylabel(r'$x_n$')
ax1.set_title('Bifurcation diagram')

ax2.plot(r, lyapunov)
ax2.axhline(0.0, color='0.4', linewidth=1)
ax2.set_xlabel(r'$r$')
ax2.set_ylabel(r'$\lambda$')
ax2.set_title('Lyapunov exponent')
ax2.grid()
fig.tight_layout()
plt.show()

The map is deterministic for every value of $r$; what changes is the character of its long-time behavior. For $0<r<1$, the trajectory approaches $x=0$. For $1<r<3$, it approaches a nonzero fixed point. At $r=3$, that fixed point loses stability and a period-doubling cascade begins. The stable period-2 orbit persists until $r=1+\sqrt{6}\approx3.4494897$, followed by a period-4 orbit until approximately $r=3.54409$, and then by further period doublings.

The period doublings accumulate at $r_\infty\approx3.56994567$. At this onset of chaos, the Lyapunov exponent approaches zero; it does not suddenly become positive at every nearby value. Beyond this point, broad chaotic regions with $\lambda>0$ are interrupted by periodic windows with $\lambda<0$. The prominent period-3 window near $r\approx3.83$ is one example.

Homework 1

For the quadratic map used to construct the Mandelbrot set, compute a finite-time orbit-stability indicator for the bounded orbits and plot it in the complex $c$ plane. This quantity is often called a finite-time Lyapunov exponent, and is closely related to Lyapunov exponent, but it is associated with an individual orbit; it is not a single Lyapunov exponent of the Mandelbrot set as a whole.

Mark points whose orbits escape separately.

Background material can be found at wikibooks.org.

The Mandelbrot set is generated by the quadratic map $$z_{n+1}=f_c(z_n)=z_n^2+c,\qquad z_0=0.$$ For a small perturbation along an orbit, $$\delta z_N\approx\delta z_0\prod_{n=0}^{N-1}f_c'(z_n),\qquad f_c'(z)=2z.$$ This suggests the finite-time average $$\lambda_N(c)=\frac{1}{N}\sum_{n=1}^{N}\log|2z_n|.$$

We begin the sum at $n=1$ because the conventional Mandelbrot orbit starts at the critical point $z_0=0$, where $f_c'(z_0)=0$. Including that first derivative would make the product zero for every $c$ and would hide the subsequent stability behavior. We first discard a transient, then average over a finite number of iterates. If the orbit escapes beyond $|z|=2$, we record NaN rather than inventing a Lyapunov value. Near the boundary, the resulting picture depends on the transient and averaging lengths, which is itself an important finite-time effect.